library(plotly)
library(readr)
library(gganimate)
library(ggmap)
library(dplyr)
library(ggplot2)
library(highcharter)
library(ISOcodes)
library(countrycode)
library(stringr)
library(lubridate)
worldwide = read_csv("data_eq/worldwide.csv")
worldwide$year = format(worldwide$time, "%Y")
for (i in 1:nrow(worldwide)) {
worldwide$place[i] = worldwide$place[i] %>% str_split(", ") %>% .[[1]] %>% tail(1)
}
historical = read_rds("data_eq/historical_web_data_26112015.rds")
disaster = read_delim("data_eq/disaster.txt", "\t", escape_double = FALSE, trim_ws = TRUE)
iran_equake = read_rds("data_eq/iran_earthquake.rds")با استفاده از داده های زلزله ها در ایران و جهان به سوالات زیر پاسخ دهید.
۱. با استفاده از داده های historical_web_data_26112015.rds و استفاده از نمودار پراکنش سه بعدی بسته plotly نمودار طول، عرض و عمق زلزله ها را رسم نمایید. علاوه بر آن بزرگی هر نقطه را برابر بزرگی زمین لرزه قرار دهید.
plot_ly(historical, x = ~Longitude, y = ~Latitude, z = ~Depth, size = ~Magnitude,
marker = list( symbol = 'circle', sizemode = 'diameter'), sizes = c(1, 30))۲. پویانمایی سونامی های تاریخی را بر حسب شدت بر روی نقشه زمین رسم نمایید.(از داده زلزله های بزرگ استفاده نمایید.)
disaster %>%filter(FLAG_TSUNAMI == "Tsu" & !is.na(LATITUDE) & !is.na(LONGITUDE) & !is.na(LOCATION_NAME)& !is.na(EQ_PRIMARY)) %>% arrange(EQ_PRIMARY) -> t1
ggplot() + geom_polygon(data = map_data("world"), aes(long, lat, group=group), fill = "white", color = "lightblue") +
geom_point(data = t1, aes(x = LONGITUDE,
y = LATITUDE,
frame = YEAR,
color = EQ_PRIMARY,
size = EQ_PRIMARY)) +
scale_color_continuous(low = "yellow", high = "red", guide = F) +
scale_size(guide = F) -> p
p#gganimate(p, filename = "t.gif")gganimate مشکل داشت و خروجی نمیداد!!!!
۳. نمودار چگالی دو بعدی زلزله های تاریخی ایران را رسم کنید.( از داده iran_earthquake.rds و لایه stat_density_2d استفاده نمایید).
tehran_map = read_rds("data_eq/Tehrn_map_6.rds")
ggmap(tehran_map) + stat_density_2d(data = iran_equake , geom = "polygon", aes(x = Long , y = Lat , fill = ..level.. , alpha = ..level..)) + scale_alpha(range = c(0.4, 0.75), guide = FALSE) + scale_fill_gradient(low = "green", high = "red", guide = FALSE)۴. احتمال اینکه در ایران در پنج سال آینده زلزله به بزرگی هفت ریشتر رخ دهد را محاسبه کنید. (از احتمال شرطی استفاده کنید.)
۵. بر اساس داده های زلزله های بزرگ ابتدا تعداد و متوسط کشته زلزله ها را بر حسب کشور استخراج نمایید. سپس نمودار گرمایی تعداد کشته ها را بر روی کره زمین رسم نمایید.(مانند مثال زیر!)
disaster %>% filter(!is.na(TOTAL_DEATHS), !is.na(COUNTRY)) %>% group_by(COUNTRY) %>% summarise(TOTAL = sum(TOTAL_DEATHS), AVE = mean(TOTAL_DEATHS)) %>% mutate(COUNTRY_CODE = countrycode(COUNTRY , "country.name", "iso3c")) -> t4
hcmap(data = t4,joinBy = c("iso-a3", "COUNTRY_CODE"),value = "TOTAL",name = "TOTAL DEATHS")hcmap(data = t4,joinBy = c("iso-a3", "COUNTRY_CODE"),value = "AVE",name = "AVERAGE DEATHS")۶. با استفاده از داده لرزه های بزرگ و به وسیله طول، عرض، شدت، عمق مدلی برای پیش بینی تعداد کشته های زلزله بیابید.
glm(TOTAL_DEATHS ~ LONGITUDE + LATITUDE + FOCAL_DEPTH + EQ_PRIMARY, disaster, family = "poisson") %>% summary()##
## Call:
## glm(formula = TOTAL_DEATHS ~ LONGITUDE + LATITUDE + FOCAL_DEPTH +
## EQ_PRIMARY, family = "poisson", data = disaster)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -384.38 -61.46 -39.20 -21.18 1521.82
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.389e+00 5.064e-03 -274.3 <2e-16 ***
## LONGITUDE 8.578e-04 6.032e-06 142.2 <2e-16 ***
## LATITUDE 1.973e-02 2.745e-05 718.9 <2e-16 ***
## FOCAL_DEPTH -2.631e-02 4.480e-05 -587.3 <2e-16 ***
## EQ_PRIMARY 1.348e+00 6.850e-04 1968.1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 18100341 on 1051 degrees of freedom
## Residual deviance: 13160953 on 1047 degrees of freedom
## (4974 observations deleted due to missingness)
## AIC: 13166049
##
## Number of Fisher Scoring iterations: 8
۷. با استفاده از داده worldwide.csv به چند سوال زیر پاسخ دهید. تحقیق کنید آیا می توان از پیش لرزه، زلزله اصلی را پیش بینی کرد؟
worldwide$date = reshape2::colsplit(worldwide$time, pattern=" ", names=c("Part1", "Part2"))[,1]
worldwide %>% group_by(place, date) %>% mutate(mainmag = max(mag), ismain = (mag == max(mag))) %>% arrange(place,date,-mainmag) %>% mutate(timedif = as.numeric(first(time)-time)) %>% filter(timedif<0)->t7
sampleindex = sample(nrow(t7),nrow(t7)*0.9)
train = t7[sampleindex,]
test = t7[-sampleindex,]
fit1 = lm(data=train,mainmag~depth+mag)
predict1 = predict(fit1,test)
summary(fit1)##
## Call:
## lm(formula = mainmag ~ depth + mag, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.3577 -0.5962 -0.2049 0.4144 4.0499
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.176e+00 2.084e-02 104.41 <2e-16 ***
## depth -8.814e-04 4.976e-05 -17.71 <2e-16 ***
## mag 6.730e-01 5.592e-03 120.35 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8615 on 32536 degrees of freedom
## Multiple R-squared: 0.3085, Adjusted R-squared: 0.3085
## F-statistic: 7258 on 2 and 32536 DF, p-value: < 2.2e-16
summary(predict1-test$mainmag)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -3.78070 -0.41130 0.19510 -0.00004 0.59146 1.35565
fit2 = lm(data=train,timedif~depth+mag)
predict2 = predict(fit2,test)
summary(fit2)##
## Call:
## lm(formula = timedif ~ depth + mag, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -55821 -18562 2446 20702 37117
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -43383.176 565.162 -76.762 < 2e-16 ***
## depth -4.111 1.349 -3.047 0.00231 **
## mag 2580.180 151.633 17.016 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 23360 on 32536 degrees of freedom
## Multiple R-squared: 0.008823, Adjusted R-squared: 0.008762
## F-statistic: 144.8 on 2 and 32536 DF, p-value: < 2.2e-16
summary(predict1-test$timedif)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 6.04 13964.08 32141.07 34521.58 53215.76 85840.84
۸. گزاره " آیا شدت زلزله به عمق آن بستگی دارد" را تحقیق کنید؟ (طبیعتا از آزمون فرض باید استفاده کنید.)
worldwide %>% filter(!is.na(mag),!is.na(depth)) ->t8
cor.test(t8$depth,t8$mag)##
## Pearson's product-moment correlation
##
## data: t8$depth and t8$mag
## t = 50.03, df = 60629, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1914584 0.2067469
## sample estimates:
## cor
## 0.1991147
ggplot(t8 %>% head(1000),aes(x = depth,y = mag)) + geom_point()۹. میانگین سالانه زلزله ها را بر حسب کشور به دست آورید. آیا میتوان دلیلی در تایید یا رد تئوری هارپ ارائه کرد.
worldwide %>% group_by(place,year) %>% summarise(avemag = mean(mag,na.rm = T))## # A tibble: 925 x 3
## # Groups: place [?]
## place year avemag
## <chr> <chr> <dbl>
## 1 Afghanistan 2015 4.33
## 2 Afghanistan 2016 4.35
## 3 Afghanistan 2017 4.30
## 4 Afghanistan 2018 4.35
## 5 Alabama 2016 2.75
## 6 Alabama 2017 2.66
## 7 Alaska 2015 2.99
## 8 Alaska 2016 3.08
## 9 Alaska 2017 3.11
## 10 Alaska 2018 3.18
## # ... with 915 more rows
ggplot() + geom_polygon(data = map_data("world"), aes(long, lat, group=group), fill = "white", color = "lightblue") +
geom_point(data = worldwide %>% filter(mag > 4), aes(x = longitude,
y = latitude,
color = mag,
alpha = .1)) +
scale_color_continuous(low = "yellow", high = "red", guide = F) +
scale_size(guide = F)۱۰. سه حقیقت جالب در مورد زلزله بیابید.
آمار مرگ در اثر سونامی از زلزله بیشتر است.
wilcox.test(
disaster %>% filter(!is.na(FLAG_TSUNAMI) , !is.na(TOTAL_DEATHS)) %>% .$TOTAL_DEATHS,
disaster %>% filter(is.na(FLAG_TSUNAMI) , !is.na(TOTAL_DEATHS)) %>% .$TOTAL_DEATHS,
alternative = "greater"
)##
## Wilcoxon rank sum test with continuity correction
##
## data: disaster %>% filter(!is.na(FLAG_TSUNAMI), !is.na(TOTAL_DEATHS)) %>% and disaster %>% filter(is.na(FLAG_TSUNAMI), !is.na(TOTAL_DEATHS)) %>% .$TOTAL_DEATHS and .$TOTAL_DEATHS
## W = 317630, p-value = 8.633e-11
## alternative hypothesis: true location shift is greater than 0
مرگ و میر در ساعات کاری کمتر است.
disaster %>% filter(!is.na(HOUR)) %>% group_by(HOUR) %>% summarise(count = n(),avemag = mean(EQ_PRIMARY,na.rm = T),avedeath = mean(TOTAL_DEATHS,na.rm = T)) -> t11
sum = sum(t11$count)
ggplot() + geom_bar(data = t11,aes(x = HOUR,y = avedeath),stat = "identity")